Dynamical QCD simulations on anisotropic lattices 
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The simulation of QCD on dynamical {Nf = 2) anisotropic lattices is described. A method for 
nonperturbative renormalisation of the parameters in the anisotropic gauge and quark actions is 
presented. The precision with which this tuning can be carried out is determined in dynamical 
simulations on 8^ x 48 and 8"* x 80 lattices. 



I. INTRODUCTION 

The advantages of simulations with anisotropic lat- 
tices are well understood and the method has been 
used for precision determinations of an extensive range 
of quantities in the quenched approximation to QCD 
[E Hi Hi Hi IE ly| ' In general a 3+1 anisotropy is employed 
where the lattice spacing in the temporal direction, at, 
is made fine whilst keeping the spatial lattice spacing as 
relatively coarse. The advantages of this approach are 
two-fold. The improved resolution in the temporal direc- 
tion means that states whose signal to noise ratio falls 
rapidly can be more reliably determined. The high com- 
putational cost of this improvement is offset by savings 
in the coarse spatial directions. 

The isotropic lattice (whose spacing in the four space- 
time directions is ax = ay = az = at = a) regulates 
QCD in a way that breaks the continuous Euclidean 
symmetry down to the finite group of rotations of the 
hypercube. Luckily the relevant operators that trans- 
form trivially under these two groups are the same and 
so there is no renormalisation of the speed of light on 
the isotropic lattice. Once an explicitly anisotropic lat- 
tice action is introduced with a^ = ay ~ a^ = as and 
at 7^ fls, the rotational symmetry of the theory is the 
cubic point group. For the gluons, there are now two 
distinct operators not related by rotations at dimension 
four: {Tr E'^,Tt B^]; while for the quarks the set of 
dimension four operators {'iplp'ipi rn'ipipj grows to a set 
with three members: {^'ip'yiDiip, -iPjqDoijj, rwipip}- As a 
result, two new parameters appear in the action, and for 
the continuum limit to represent QCD these parameters 
must be determined such that a physical probe of the 
vacuum at scales well below the cut-off appears to have 
full Euclidean symmetry. The nonperturbative determi- 
nation of these extra action parameters is the subject of 
the present paper. 

In quenched QCD the anisotropy in the gauge sec- 
tor, ^g, and the quark sector, f,, can be tuned inde- 
pendently and post hoc using two separate criteria. The 
precision and mass-dependence of the determination of 
^g was investigated for the action we use in Ref. 7]. It 
was found that this parameter could be determined at 
the percent level from the energy-momentum dispersion 
relation. The mass dependence was found to be mild for 



quark masses in the range m^ < mq < rric when the tun- 
ing was carried out at the strange quark mass, ttIj. In 
Refs. la, 13, a determination of the gluonic parameter was 
made to similar precision. 

We would like to use anisotropic lattices in simulations 
with Nf = 2 for realistic phenomenologically-relevant 
calculations. In dynamical QCD the tuning procedure 
becomes more complicated because of the interplay be- 
tween the quark and gluon sectors and the parameters 
must be simultaneously determined. There are several 
issues to resolve. Firstly, can this simultaneous tuning 
be accomplished; secondly, to what precision is the renor- 
malised anisotropy determined; and thirdly, what is the 
mass-dependence of the renormalised anisotropy. Here 
we will focus on the first two issues, and leave the ques- 
tion of the mass dependence to a later study. 

The paper is organised as follows. Section IITI gives the 
details of the gauge and quark actions used in this investi- 
gation. Section llTTl describes the tuning methodology and 
is followed in Section HVI bv the results for the values of 
the tuned bare (input) parameters ^" and ^°. Section [V] 
contains our conclusions and future plans. 



II. THE ACTION AND PARAMETERS 

We begin with a brief description of the anisotropic 
action used in this study. The details of the tuning pro- 
cedure described in Section IIIII do not depend on the 
specific action used. Further description of the action 
can be found in [7| where the tuning for the same action 
in the quenched approximation was discussed. 

The gauge action is a two-plaquette Symanzik- 
improved action [l3| previously developed for high- 
precision glueball studies and given by 
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where ilg and fit are spatial and temporal plaquettes. 
ilf and fl^ are 2x1 rectangles in the (i,j) and (i,t) 
planes respectively, il^* is constructed from two spatial 
plaquettes separated by a single temporal link. Us and 
Ut are the mean spatial and temporal gauge link values 



respectively. The action has leading discretisation errors 
of 0{al,at,asas). 

For fermions an action specifically designed for large 
anisotropics is used. The usual Wilson term removes 
doublers in the temporal direction whereas spatial dou- 
blers are removed by the addition of a Haniber-Wu term. 
The action has been described in detail in Ref. 7] and 
has leading classical discretisation errors of 0{atraq). In 
terms of continuum operators, it can be written 
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which highlights the different treatment of temporal and 
spatial directions, r is the usual Wilson coefficient which 
is applied in the temporal direction only in this action 
and is set to unity. The analagous parameter in the spa- 
tial directions is s, which parameterises a term that is 
irrelevant in the continuum limit. A precise tuning of 
this parameter is not necessary: in practice we choose 
s = 1/8, so that the energy of a propagating quark at tree 
level increases monotonically across the Brillouin zone. 
Stout-link smearing [ij was used for the gauge fields 
in the fermion matrix. Two stoutening iterations were 
used, with a parameter p = 0.22. This was fixed for all 
simulations, and chosen to approximately maximise the 
expectation value of the spatial plaquette on the stout 
links. 

This study was carried out on ^ x 48 and 8^ x 80 
anisotropic lattices with a spatial lattice spacing a^ « 
0.2fm and a target anisotropy ■^ = 6. The bare sea quark 
mass was set to atra^ — —0.057 in all runs. A set of 
gauge configurations, distributed across ten independent 
Markov chains, was generated for each set of input pa- 
rameters (Cg,Cq)- Valence quark propagators were gener- 
ated with the same mass as the sea quarks. 

To determine the statistical uncertainties, 1000 boot- 
strapped sets of configurations were taken and analysis 
was done on these bootstrapped sets. Both point and all- 
to-all propagators were used. Some preliminary results 
using point propagators on 8^ x 48 lattices were presented 
in Ref. [H. 



III. METHODOLOGY 

The bare parameters, ^„ and ^J?, are renormalised by 
demanding that physical probes exhibit euclidean sym- 
metry. In principle, any physical quantity can be used; 
however, it should be easily determined to high precision. 
In this study we have used the sideways potential and the 
pion energy-momentum dispersion relation for the gauge 
and fermion sectors respectively. 

The gauge anisotropy ^g is determined from the in- 
terquark potential i^iQj . The static source propagation is 
chosen to be along a coarse direction allowing the sources 
to be separated along both course and fine axes. The po- 
tential is determined at the same physical distance for 
these two cases. The input anisotropy is constrained so 



that the two calculations yield the same value of the po- 
tential, Vs{x) = Vt(i/^) for a target anisotropy f. For a 
given input anisotropy ^ and target anisotropy ^ we can 
determine the mismatch parameter Cg = Vs(a^)/Vt(i/^). 
If X is in the regime where the potential is nearly linear, 
the mismatch parameter is approximately related to the 
actual gauge anisotropy, Cg w fg/'C- 

The quark anisotropy can be determined from the 
pseudoscalar dispersion relation. The anisotropy is in- 
versely proportional to the square root of the slope 
of the dispersion relation and demanding a relativis- 
tic energy-momentum relation imposes a renormalisa- 
tion condition on the bare parameter ^ . The ground 
state energy i?o was determined for a range of momenta, 

22ii and we aver- 



e {0,1,2,3,4,5,6}, where p„ 



age over equivalent momentum values. The two-point 
correlator data were modelled with single exponentials 
and a ^^-minimisation was used to determine the best- 
fit ground state. These values were used to generate an 
energy-momentum dispersion relation. 

In the quenched approximation this procedure is rela- 
tively easy since ^ and ^„ can be determined indepen- 
dently. For dynamical simulations it is no longer possible 
to simply fix ^^ and then tune ^^ to a consistent value, 
since changing ^ will affect the measurement of ^g . Ex- 
plicitly, changing the value of i^ necessitates a regener- 
ation of the background fields with the new value of i^ 
which in turn will change the measured anisotropy ^g of 
the background fields. The solution to this problem is a 
simultaneous two-dimensional tuning procedure [13 . 

A linear dependence on the parameters ^g and ^g was 
assumed for a small region. Three initial sets of configu- 
rations were generated and the renormalised anisotropy 
was determined. Planes were defined for both output val- 
ues of ^g and ^q i.e. values a,/?, 7 were found to satisfy 
Ca — ota^g + /3aCg + la for the renormaliscd anisotropy 
^a,a = ffi 9 measured for each input {&,^g)- The inter- 
section of these planes with the required (target) output 
value yields the tuned point. The statistical uncertain- 
ties were determined using bootstrap resampling, with a 
common bootstrap ensemble used for all measurements. 
When more than three simulation points were available 
a plane was defined using a constrained-x^ fit. 

All observables were estimated using the Monte Carlo 
method. An ensemble of 250 gauge field configurations 
divided across 10 Markov chains was generated using the 
Hybrid Monte Carlo (HMC) algorithm |Q|. Approxi- 
mately 5000 CPU hours were needed in order to generate 
each set of configurations. The HMC algorithm can be 
used for these simulations without modification. One ob- 
servation serves to improve performance, however. HMC 
adds a set of momentum variables conjugate to the gauge 
fields, but each conjugate momentum can be added with 
a different gaussian variance without changing the valid- 
ity of the method. In isotropic simulations this is not a 
useful property, and all momentum co-ordinates are cho- 
sen to have unit variance. For the anisotropic lattice, the 
temporal and spatial gauge fields have different interac- 
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3 4 5 
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1.51 


1.528 


1.514 1.544 1.522 


e. 


6.0 


7.5 


7.5 8.72 8.83 


e. 


8.0 


7.0 


8.0 6.65 7.44 



TABLE I: Input parameters for the five dynamical simulations 
performed in this tuning procedure. The bare quark mass is 



atm„ 



-0.057 for all runs. 



tions, and different momenta become useful. If the HMC 
hamiltonian is 



H = 



1 



X^l. 
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an extra tunable parameter, fit (the variance of the 
temporal link momenta), has been added to the algo- 
rithm which can be used to optimise acceptance by the 
Metropolis test. This is equivalent to using two distinct 
integration step-sizes for the spatial and temporal degrees 
of freedom. Some brief numerical experiments suggest 
that a temporal leap-frog step-size smaller by a factor ^ 
is close to optimal, and this is borne out by considera- 
tions of free field theory. 



IV. RESULTS 



Run 



Cg = Vs{x)/Vt{t/^) at different (T,R) 
(1,3) (1,4) (2,3) (2,4) (3,3) (3,4) 



0.972(2) 0.959(3) 0.972(7) 0.965(13) 0.991(25) 1.13(8) 
0.951(2) 0.941(4) 0.945(8) 0.926(18) 0.942(34) 0.89(9) 
0.994(2) 0.990(3) 0.991(7) 0.998(13) 0.965(25) 1.01(7) 



TABLE 11: The gluon anisotropy parameter Cg for different 
separations, R and times, T. The final results were deter- 
mined from data at T = 2 and R = 3. 



B. Dispersion relations 

Pseudoscalar meson correlators were computed using 
traditional point propagators as well as all-to-all propa- 
gators [l5| with time and colour dilution and no eigen- 
vectors. 

To determine optimal fit ranges for exponential fits 
to the correlator data, sliding window (imin) plots were 
used: the correlation function was fitted in a range from 
tmin to iniax whcrc imax was fixcd to the largest value 
compatible with a good fit, and imin was varied. An ex- 
ample of such a plot is given in Fig. ^ The fit range was 
chosen so the fit would be stable with respect to small 
variations in imin- The same fit ranges and smearing pa- 
rameters were chosen for all simulation points in order 
to obtain a consistent determination of the dispersion re- 
lation. The final fit ranges are given in Table IIIII In 



The input anisotropy parameters used are given in Ta- 
ble ^ We started by choosing three points (Runs 1-3) 
in the (^° ^°) plane, and generated configurations at two 
further points as a result of the tuning procedure. The 
final tuning was performed on 8'^ x 80 lattices, using data 
from runs 1, 4 and 5 as these spanned the largest area of 
the plane. 



A. Interquark Potential 

The gluon anisotropy is determined from the static 
quark potential at a selected distance R. In practice this 
is done by determining the effective energy for the static 
quark-antiquark configuration at separation R at some 
time T. It is then important to choose values for R and 
T where the potential is well determined and the value 
obtained for Cg is stable with respect to small variations 
in R and T. The same values for R and T must then be 
used for all runs in order to have a consistent procedure. 

Table im shows Cg for different R and T, on the 8^ x 80 
lattices. We see that the values are generally quite con- 
sistent for each run. Looking more closely at the effective 
potential for each i? as a function of T, we find that it has 
not yet reached a plateau at T = 1, while the value for 
T = 3 is consistent within errors with that for T = 2. We 
choose (r, R) = (2, 3) as our optimal parameters, since 
this yields reasonably small statistical errors, while R is 
large enough to be in the linear regime. 
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FIG. 1: A typical tmin plot, showing the energy for momentum 
n^ = 1 on run 1, 8'^ x 80 lattices from fits to time ranges 
imax = 40 for various tmin- A stable ground state energy 
determination, with a good x^, is achieved for 22 < tmin < 30. 

our initial analysis data from a 8"^ x 48 lattice were used. 
However, a reliable extraction of the ground state energy 
proved difficult. In particular, it was observed that the 
energy either did not reach a plateau until near the end 
of the lattice or did not plateau at all. To resolve this 
problem the simulation was repeated on a longer, 8^ x 80 
lattice. An immediate improvement in the quality of the 
fits was observed. The ground state energy was deter- 
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40 
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24 


40 
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21 


40 
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19 


40 



TABLE III: Fit ranges. 



mined from fits over at least 15 timeslices and was stable 
with respect to changes in imin- The effect of the longer 
lattice is illustrated in Figure |21 This plot also com- 
pares simulations using point and all-to-all propagators. 
The all-to-all propagators lead to improved precision in 
the fitted energies. The central values are in agreement 
with the energies determined using point propagators but 
the statistical error is smaller. The final tuned parame- 



8'x48 point 
8'x80 point 
8'x80 all-to-all 



a~ 0.04 



0.02 





FIG. 3: Dispersion relations from runs 1, 4 and 5 on 8 x 80 
lattices using all-to-all propagators. The solid line is a fit 
to the four points and the dotted fines are the 68% confi- 
dence levels. The quality of all three fits is very good with 
X^/Nd.f. = 2.0/2, 1.9/2, 2.0/2 for runs 1,4 and 5 respectively. 





8^ x48 


8^ X 80 


Run 


•^9 ?9 


Cg ?q 


1 
2 
3 
4 
5 


0.991(3) 4.98(6) 
0.986(3) 6.27(4) 
1.001(3) 5.18(6) 
0.985(5) 6.47(5) 
0.995(3) 5.80(5) 


0.972(7) 5.54(6) 

0.945(8) 7.08(5) 
0.991(7) 6.95(8) 



TABLE IV: Table of measured output anisotropics at each of 
the run points. The errors are statistical only. 



FIG. 2: A comparison of the dispersion relations determined 
from an 8^ x 48 lattice and an 8^ x 80 lattice. The solid 
lines are the best fits and the dotted lines are the 68% confi- 
dence levels. The figure also shows a comparison of all-to-all 
propagators and point propagators on the same (longer) lat- 
tice. The plot shows that the ground state energies have not 
reached a plateau on the shorter lattice. On the longer lattice 
the all-to-all and point data agree, while higher precision is 
achieved with all-to-all propagators. 

ters were determined using all-to-all propagators on the 
8'^ X 80 lattices. We find consistently good fits for all 
runs for the first four momenta considered {n? = 0, 1,2 
and 3). The renormalised quark anisotropy is therefore 
determined from fits to these momenta. Figure 13 shows 
the pseudoscalar dispersion relations for Runs 1, 4 and 5 
which are used to determine the tuned point. 



C. Plane fits 

Table HVl shows the output anisotropics determined on 
the 8*^ X 48 and 8'^ x 80 lattices for the five simulation 
points. As a check on the stability of our tuning proce- 
dure, we have repeated the calculation using different 
values of R and T in the determination of the gluon 



anisotropy. The results are shown in Fig. ^ The plot 
shows that the anisotropics are insensitive to a change in 
R but that increasing the value of T from two to three 
leads to large statistical uncertainty, particularly in the 
gluon anistropy. For these reasons we choose R — 3 and 
T = 2 for our analysis. 



D. Simulation with tuned parameters 

Applying the plane fit procedure of Sec. lIVCl to a sub- 
set of configurations of Runs 1, 4 and 5 we obtained pre- 
liminary, tuned parameters ^^ — 8.O6I7, ^° = 7.52I15. 
250 configurations were generated with these parameters, 
and Cg and ^^ determined using the same values for R, 
T and fit ranges as in Sections IIV Al and IIVBI We find 
Cg = 0.983(6), Cg = 6.21(9). We see that both quark and 
gluon anisotropics are within 3% of the target value of 6. 
Although the anisotropics are not equal within statistical 
errors, we note that there are still systematic uncertain- 
ties at the percent level, in particular for ^g, as shown in 
Table Hll For example, if we choose i? = 3, T = 3 we find 
c, = 1.01(2). 

We repeated the plane fit procedure including the new 
information from Run 6. Figure [S] shows the resulting 
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FIG. 4: Tuned values of input parameters (^j, ^q) determined 
from the plane fit procedure on the 8^ x 80 lattice. The plots 
show the results for different values of R and T used to de- 
termine the gluon anisotropy. Each point corresponds to one 
bootstrap sample. 
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FIG. 6: The potential between fundamental static color 
sources for run 6, measured from static propagation in a 
coarse direction. Lines show fits to the Cornell potential, and 
are used in a crude determination of the lattice spacing. 



V. CONCLUSIONS 
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FIG. 5: As in Fig 0] The figure shows the results from a 
plane fit using parameters from runs 1, 4, 5 and 6 (marked 
with an x). The big red (gray) cross at {i°,i°) = (8.42, 7.43) 
indicates the result of the best fit. 



scatterplot determined on the 8^ x 80 lattice from runs 
1, 4, 5 and 6. The intersection points shift in a direction 
to move Cg and S,q even closer to the target anisotropy. 
In order to get a rough idea of the physical scales of 
these lattices, we compute the pion mass, the rho mass 
and the string tension. We find atw^ = 0.066(1) and 
atnip = 0.120(5), which gives rriT^/mp = 0.54, while a 
crude measurement (shown in Fig.|SJ) of the string tension 
gives as = 0.18fm. A more precise determination of the 
lattice spacing will be obtained from the IP-IS splitting 
in charmonium [l6j. 



We have performed a first simulation of 2-flavour QCD 
with improved Wilson fermions on anisotropic lattices, 
with both quark and gluon anisotropics tuned to ^ = 6 
[23 . The tuning was based on a linear Ansatz for the de- 
pendence of renormalised anisotropics on bare anisotropy 
parameters in a region of parameter space. The results 
from the final run demonstrate that the tuning proce- 
dure, described in Sec. IIIII works satisfactorily. 

The final, tuned point was found to lie marginally out- 
side the triangle used for the plane fit procedure, so the 
end result was based on an extrapolation rather than 
an interpolation. This increases both the statistical and 
systematic uncertainties of the determination. To avoid 
this problem, it is important to choose a large enough 
triangle to start with, so that successive parameter de- 
terminations are always based on interpolations. 

We also found that the original (8'^ x 48) lattices used 
were too short in the time direction to allow a reliable de- 
termination of ground state energies, which were found 
to be systematically high, in particular for higher mo- 
menta. This led in turn to systematically high values for 
^q. The adoption of lattices with longer time extent was 
a crucial step in the procedure. As Table UTTl shows, the 
optimal fit ranges were generally found to be beyond the 
range of the shorter lattice. 

We were able to determine the tuned parameters 
(^°,^") with a statistical uncertainty of 1% and 3% re- 
spectively from our ensembles of 250 configurations. In 
addition, there are three main sources of systematic un- 
certainties: 

1. The R and T values used in the determination of 
the sideways potential, and the fit ranges used in 
the determination of the pseudoscalar dispersion re- 
lation. Since the fit ranges are chosen to give stable 



ground state energies, we can safely assume that 
the latter is a small effect. The effect of varying 
R is also small, as shown in Fig. ^ There may 
be a systematic error arising from the choice of T, 
but this is obscured by the larger statistical uncer- 
tainties in the T = 3 data, particularly in the ^^ 
direction. 

2. Lattice sizes. The pion dispersion relation is un- 
likely to be strongly affected by the finite lattice 
volume, but the static quark potential may contain 
finite volume errors which affect our results. We 
will be performing simulations at the tuned point 
on larger volumes, which will show whether this is 
a significant issue. 

3. Nonlinearities in the dependence of {£,g,£,q) on 
(^gjCg)- Our final fit to four points shows no ev- 
idence of any significant nonlinearity. If this were 
found to be a serious issue in any future simula- 
tion, a two-step procedure may be adopted where 
a smaller triangle centred on the preliminary tuned 
point is used in the second step. 

We have yet to verify that we get the same quark 
anisotropy from other hadronic probes, for example the 
vector meson. Differences in the anisotropies can arise 



from lattice artefacts and can thus be considered part of 
the finite lattice spacing errors. 

These lattices will in the future be employed for a wide 
range of physics investigations, including charm physics 
and heavy exotics 16], spectral functions at h igh temper- 
ature |l3, static-light mesons and baryons [Ig, strong 
decays and flavour singlets including glueballs. These 
studies will be carried out on larger lattice volumes. Sim- 
ulations on finer lattices will necessitate a new nonper- 
turbative tuning process like the one performed here; this 
will be desirable in the longer term. 
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